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Abstract 



It is well-known from the work of ISchonbucheij ()2005() that the marginal laws of a 



loss process can be matched by a unit increasing time inhomogeneous Markov process, 
whose deterministic jump intensity is called local intensity. The Sto c hastic Local Intensity 



(SLI) models such as the one proposed bv lArnsdorf and Haloerinl ()2008l) allow to get a 



stochastic jump intensity while keeping the same marginal laws. These models involve a 
non-linear SDE with jumps. The first contribution of this paper is to prove the existence 
and uniqueness of such processes. This is made by means of an interacting particle system, 
whose convergence rate towards the non-linear SDE is analyzed. Second, this approach 
provides a powerful way to compute pathwise expectations with the SLI model: we show 
that the computational cost is roughly the same as a crude Monte-Carlo algorithm for 
standard SDKs. 
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1 Introduction 



In equity modelin g;, a ma j or con cern is to get a model that fits option data. It is well-known 
from the work of iDupire that basically European options can be exactly calibrated 

by using local volatility models cr{t,x). However, local volatility models are known to have 
some inadequacy to describe real market s. To get richer dynamics , Stoch ast ic Local V o latilit y 
(SLV) models have been introduced (see Alexander and Nogueira ( 20041 ) or Piterbarg ( 20061 )) 
and consider the following dynamics for the stock under a risk-neutral probability: 

dSt = rStdt + /{YtHt, St)StdWt, 

where is an adapted stochastic process. Typically, {Yt,t > 0) is assumed to solve an 
autonomou s one din i ension al SDE whose Brownian motion may be correlated with W. From 
the work of Gvongv ( 19861 ). we know that under mild assumptions, the following choice 



r]{t,x) 



a{t, x) 



y'E[f{Yt)^\St=x] 



ensures that St has the same marginal laws as the local volatility model with a{t,x), which 
automatically gives the calibration to European option prices. This leads to the following 
non-linear SDE 

dSt = rStdt + —=J^^=a(t, SAStdWt. 

Here, we stress that the law of (If, St) steps into the diff u sion t erm. Unless for trivial choices 
of Yt and despite some attempts ( Abergel and Tachet ( 20ld )). getting the existence and 
uniqueness of solutions for this kind of SDE remains an open problem. Also, from a numerical 
perspective, the simulation of SLV models is not easy, precisely because of the computation 
of the conditional expectation. 

In this paper, we propose to tackle a very analogous problem arising in credit risk modeling. In 
all the paper, we will work under a risk-neutral probabilistic filtered space (0, {J-t)t>o, J',^)- 
As usual, J-t is the c-field describing all the events that can occur before time t and J- 
describes all the events. We consider M S N* defaultable entities (for example, M = 125 for 
the iTraxx). We assume that all the recovery rates are deterministic and equal to 1 — Lgd 
for all the firms within the basket. The loss process {Lt,t > 0) is given by 



Lt 



Lgd 
M 



Xt, > 



where Xt is the number of defaults up to time t. Clearly, X takes values in Cm '■= {0, . . . , M}. 
Thanks to the assumption of deterministic recovery rates, and under the assumption of de- 
terministic short interest rates, it is well known that CDO tranche p rices only depe nd on the 
marginal laws of the loss process (see for example Remark 3.2.1 in Alfonsi ( 201ll )). Let us 



assume then for a while that we have found marginal laws {P{Xt = k),k € ^M)te[o,T] which 
perfectly fit C DO tranche p r ices u p to maturity T > 0. Then, under some mild assumptions, 
we know from ISchonbucheri (j2005l ) that there exists a non- homogeneous Markov chain with 
only unit increments which exactly matches these marginal laws. Somehow, this result plays 
the same role for the loss as Dupire's result for the stock. 

The loss model obtained with a non-homogeneous unit-increasing Markov chain is known 
in the literature as local intensity model. It is fully described by the local intensity A : 
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M+ X Cm which gives the instantaneous rate of having one more default in the basket. 

In the sequel, we will assume that the local intensity A : M+ x Cm 1^+ has been calibrated 
to market data and perfectly matches Index and CDO tranche prices. We make the following 
assumptions: 

• Mx £ Cm-.^ A(t, x) is a cadlag function, 

• > 0, A(t,M) = 0. 

In particular, we have A = mayix^cM sup^gp.T] -^(^) ^) < ^xd. In this setting, the instantaneous 
jump rate at time t from x to x + 1 is given by \{t—,x). Thus, the local intensity model 
corresponds to a time-inhomogeneous Markov chain making unit jumps with this rate. 
One may like however get richer dynamics than the ones given by the local intensity model. 
Then, we can proceed in the same way as for the stochastic volatility models. Let us consider 
{Xt)t>o a general (7"t)-adapted cadlag real process, a function / : R — )■ M+ and a function 
rj : ]R_|_ X Cm ^ ^+ satisfying the same assumptions as A. We assume that the default 
counting process {Xt,t > 0) has jumps of size 1 with the rate 

r,it-,Xt.)f{Yt^). 

By analogy with the equity, w e name this kind of mo del a Stochastic Local Intensity (SLI) 
model. Then, it is known (see Cont and Minca ( 20131 )) that the local intensity model with 



the Local Intensity (LI) 7y(t— , x)E[/(5^_)|X(_ = x] has the same marginal laws as Xf. Thus, 
the SLI model will be automatically calibrated to CDO tranche prices if one takes: 

Vt > 0,x G Cm, ?7(t-,x) - Kt-,x) 



This approach has been used in the literature by 



E[f{Yt_)\Xt-=xy 
Arnsdorf and Halperin ( 20081 ). and 



slightly different way bv iLopatin and Misirpashaevi (|2008l ). However, up to our knowledge 



m a 



there is no proof in the literature of the existence nor uniqueness of such a dynamics. 

The first scope of this paper is to solve this problem. At this stage, we need to make our 

framework precise. We assume through the paper that: 

/ : M R is continuous, s.t. Vx G M, < / < /(x) < / < oo. (1.1) 

We assume that the probability space {ft, J-, P) contains a standard Brownian motion (Wt, t > 
0), a sequence of independent uniform random variables {U^)keN, and a sequence (ii'")„gN 
of independent exponential random variables with parameter We set = Yl^=i^^ 

for A: G N*. The random variables [T^ ,U^)k&^* will enable us to define a non homogeneous 
Poisson point process with jump intensity r/(i— , Xf_)/(1^„). We are interested in studying 
the following two problems in which we assume that Yj is a process with values either in N 
(discrete case) or in M (continuous case). In the discrete case, we are interested in finding a 
predictable process {Xt.,Yt)t>Q such that 



-Tf E|/CS'j,fc_)|Xyfe_ 



^0 = yO) and for each A; > 0, (l(,t G [Tfc,Tfc+i)) is a continuous time Markov chain 

Xt _ ^1 
'ij - f^ij 



with transition rate fif/ = ji- 



(1.2) 
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For the sake of simplicity, we consider in the discrete setting that X and Y do not jump 
together almost surely. 

In the continuous case, the corresponding problem is to solve the following stochastic differ- 
ential equation: 



Xt - Xq + Ylk,Tk<t 1 



(1.3) 



\Yt = y^ + £b{s,Xs,Ys)ds + £a{s,Xs,Ys)dWs + £ ^{s-, Xs-,Ys^)dXs, 



for xq G £a/> yo £ ^ and given real functions b, a and 7. Th is frai nework embeds in particu- 
l ar th e dynamics suggested by Arnsdorf and Halperin ( 20081 ) and Lopatin and Misirpashaev 
( 20081 ). Under some rather mild hypotheses on Hij, b, a and 7, which will be specified in 
the corresponding sections, we will show that the above two equations admit a unique solu- 
tion. In the discrete case, we are able to show that the corresponding Fokker-Planck equation 
has a unique solution. This can be achieved by writing the Fokker-Planck equation as an 
ODE which can be studied directly. In the continuous case, this approach can hardly been 
extended: the Fokker-Planck equation leads to a non trivial PDE. Instead, we solve this prob- 
lem by introducing an interacting parti cle system. This tec hnique is know n to be powerful 
for this type of non linear problems (see Sznitman ( 1991 ) or Meleard ( 19961 )). 
The second scope of this paper is to provide a way to compute prices under SLI models. 
Indeed, interacting particle systems are not only theoretical tools to prove existence and 
uniqueness results for such equations. They give a very smart way to simulate these 
processes, theref ore enabling us to run Mont e -Carl o algorithms. This approach has been 
recently used bv lGuvon and Henrv-Laborderd (|201ll ) for Stochastic Local Volatility models. 
For the Stochastic Local Intensity models considered in this work, the conditional 
expectation is much simpler to handle. This enables us to get theoretical results on the 
convergence and also simplifies the implementation. In fact, we show in our case under some 
assumptions that the rate of convergence to estimate expectations is in 0{1/N'^) for any 
a < 1/2, where N is the number of particles. On our numerical experiments, we even 
observe on several examples a convergence which is similar to the one of the Central Limit 
Theorem, which is rather usual for Interacting Particle Systems. Besides, we show that we 
can simulate the interacting particle system with a computational cost in 0{DN), where D 
is the number of time steps for the discretization of the SDE on Y. Thus, the computational 
cost is roughly the same as a crude Monte-Carlo algorithm for standard SDEs with N 
samples. 



The paper is organized as follows. First, we study the case where Y has discrete values; this 
framework enables us to settle the problem and solve it by rather elementary tools. This part 
is independent from the rest of the paper. Second, by means of a particle system approach, 
we investigate the case where Y is real valued jump diffusion. Finally, we carry out numerical 
simulations highlighting the relevance of the particle system technique to compute pathwise 
expectation of the Process (|1.3p . 



2 The SLI model when Y takes discrete values 

The goal of this section is to prove the existence of a process {Xt, Yt)t>o satisfying (|1.2p . Unlike 
the continuous case (jl.Sp . we can get this result by elementary means, without resorting to 
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an interacting particle system. To do so, we write the Fokker-Planck equation associated to 
the process {Xt,Yt)t>o, which should be satisfied by F{Xt = i,Yt = j) for G Cm x N. 
We have 



dtp{t,i,j) = Efc^j l^kjPit, i, k) + l{i>i} f(3)p{t, i - 1, j) 

P(.0,i,j) = V(i,i) / (xo,?/o), 



where p is a function from y. Cm x N to M and 



Yl%of{j)p{t^hj) 

TX^oPit^i.j) 



If we manage to prove that the Fokker-Planck equation admits a unique solution p such that 

(2.1) 

(2.2) 



V«,jE£A/xN,Vi>0, p{t,i,j) >0, 

M oo 

Vt>0 5^J^p(t,i,j) = l, 



j=0 j=0 

then we will get that the law of a process (Xt,Yj)t>o satisfying (|1.2p is unique. Besides, we 
will also get the existence of such a process. It is easy to check that a continuous Markov 
chain {Xt,Yt)t>Q starting from (2:0,1/0) with transition rate matrix 

/^(jiji),(j2,i2)(*) = lji=j2 (li2=ji+i - li2=n) '^^,'^^1, V/'* + In=i2/^}lj2' ^i'^2 G CM,ji,j2 e N, 
where p is the solution of {£) satisfies (|1.2p . In fact, the Fokker-Planck equation of this process 



dtq(t,ij) 
q{0,xo,yo) 



Ek^jfJ-ijlit^hk) + l{i>i}^gi^/(i)g(i,« - 1, j) 
V(i,j) / (2;o,yo), 



1. 



is linear and clearly solved by p, which gives q = p. 
2.1 Assumptions and notations 

In this part, we assume that the transition rates satisfy the following hypothesis. 
Hypothesis 1 The intensity matrices {fJ'ij)i,j>o for k E Cm satisfy the following conditions: 

• V/c G Cm, Vi, j G N X N such that i / j fi'^j > 0, 

• yk£ Cm, Vi G N < 0, 

. yk G Cm, Vi G N Er=o4 = (then T.T=o,,^^^4, = -f-^u)- 
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Moreover, we assume \/k G Cm, supjgp^ I/j-^J < oo. 

We also introduce specific notations used in the discrete case. 

• We define tlie set E as the set of real sequences indexed by i G Cm and j G N: 

E ■.= {u = (n})o<i<A/,o<i : u] G M}, 
E+ := {u = (u*)o<i<M,o<i : W- G M+} and E^ := {u = (ttj)o<i<M,o<i : u) G 1^+}- 



For u e E,we set |u| = J2iio l^jl 



2.2 Solving the Fokker- Planck equation 

To be more concise, we rewrite the Fokker-Planck Equation {£) and Conditions ()2.ip — (|2.2p 
by using a sequence of functions. Let P := (-fj )o<j<Af j>o denote a sequence such that each 
Pj is a function from M+ to M. Solving {£) under Constraints (|2.ip — (|2.2p boils down to 
solving 

=^{t,P{t)), 
P{0) = Po 

under the constraints \/t > 0, P{t) > and \P{t)\ = 1. The sequence Pq is such that {Po)] = 
for 7^ (xo, yo) and (fo)^" = 1- ^ is an application from R+ x — )• ii^ given by 



(^0 



k>0 



where := ^j^^'^^'^'^ 

Remark 1. ^ zs defined without ambiguity on E^. When x G S+, difficulties may arise 

z\ 

when for some fixed i, x\ = Vj. In this case, we still have lim — = smce 

zS-B* ,2-5.a: if[Z, i) 

f ^ 93(2, i) < f for z G Thus, we can extend ^ by continuity on E^. 

We aim at solving (£) in the set of summable sequences (compatible with Condition ()2.2p ) 
and get the following result. 

Theorem 2. Equation {£') admits a unique solution on M+ satisfying Vt > 0, P{t) > and 
\P{t)\ = 1. 

To do so, we first focus on the following differential equations: 

'P'{t) =^{t,{P{t))+), 
P(0) =Po. 

Proposition 3. Equation {£") admits a unique solution on Moreover, the solution 

satisfies Vt > 0, P{t) > and \P{t)\ = 1. 

The proof of this proposition consists in first showing the Lipschitz property which gives the 
existence and uniqueness of P and then proving that Vt > 0, P{t) > and |-P(i)| = 1. This 
proof is postponed to Appendix lAl 

Then, the proof of Theorem [2] becomes obvious. The unique solution P of {£") clearly solves 
{£'), and any solution Q of {£') such that Vi > 0, Q{t) > and = 1 also solves {£") 

and thus coincides with P. 
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3 The SLI model when Y is real valued. 

3.1 Setting and main results 

We are interested in proving the existence of a process {Xt,Yt)t solving the stochastic dif- 
ferential equation p.3|) . More precisely, we will consider the following stochastic differential 
equation, 



Xt- Xo + Y.kTf'Kt 1 



^ -a7 E[/(iyfc_)|Xyfc_l [ (3.1) 



Yt = YQ + /o* 6(s, X„ Ys)ds + /o a{s, X,,Ys)dWs + /q 7(s-, y,_)dX„ 

with (possibly) random initial condition [Xq^Yq) such that E[|yo|'"] < oo for any m G N. We 
will denote in the sequel Cinit the probability law of {Xq^Yq) under P. To get existence and 
uniqueness results for (j3.ip . we will make the following assumption on the coefficients. 

Hypothesis 2 1. The functions b, a, ^ : M+ x C]\j x R — > M are measurable, with sub- 
linear growth with respect to y: 

VT > 0,3Ct > 0, Vt G [0,r],x G CM,y G R,\b{t,x,y)\ + \a{t,x,y)\+\j{t,x,y)\ < CT(l+|y|). 

2. The functions b{t, x, y) and a{t, x, y) are such that for any x £ Cm, Ho ^ there exists 
a unique strong solution for the SDE 

Yt = yG+ I b{s,x,Y,)ds+ [ a{s,x,Ys)dWs, t>0. 
Jo Jo 

This property holds if we assume for example that: 

^\b{t,x,y)-b{t,x,y')\ < CT\y-y'\ 



VT > 0, 3Ct > 0, Vt G [0, T], X G Cm, 
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{t,x,y) - cr{t,x,y')\ < CT^J\y - y'\ 



3. For any x G Cm, {t,y) ^ ^{t,x,y) is cddldg with respect to t and continuous with 
respect to y, i.e. 'y{t,x,y)=lims>t,s^t,z^ylis,x,z) andlims<:t,s^t,z^ylis,x,z) exists 
ans is denoted by 'y{t—,x,y). 

To prove the strong existence and uniqueness of a process {X,Y) solving (j3.ip . we will first 
need to prove a weak existence and uniqueness result. To do so, we introduce the Martin- 
gale Problem associated with (|3.ip . We denote by P([0,T],R) the set of cadlag real valued 
functions and consider: 

E = {{x{t),y{t))i^[QX], s-i- y £ ^([0,7"],]^) and x is a nondecr. cadlag function with values in Cm}- 

This path space is endowed with the usual Skorokhod topology for cadlag processes, and with 
the associated Borelian ir-algebra. We denote by V{E) the set of probability measures on E. 
We are looking for a probability measure Q G V{E) such that Cinit is the probability law of 
[Xq.Yq) under Q and, for any (p G 



Jo ^ '^q[J\'^u)\^u\ 

+b{u,Xu,Yu)dy(t>{Xu,Yu) + ]^a^{u,Xu,Yu)dl(l,{Xu,Yu)}du (3.2) 
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is a martingale with respect to the filtration Tt = cr((X„,l^),n < t) satisfying the usual 
conditions. Here, Xt and Yj stand for the coordinate applications: 

Xf. E Cm and Yt: E -^R 

{x,y)^x{t) {x,y)^y{t), 

and Eq denotes the expectation under the probability measure Q G V{E). Similarly, for any 
Q G V{E), we denote by Eq the expectation under the probability measure Q while E 
simply denotes the expectation under the original probability measure P. 



The following Theorem is the main result of the paper. 

Theorem 4. We assume that Hypothesis holds and consider a J^Q-measurable initial con- 
dition {Xq,Yq) such that Vm G N,E[|yo|™] < oo. Then, there exists a unique probability 
measure Q G V{E) solving the Martingale Problem ()3.2p . Besides, there exists a unique 
strong solution {Xt,Yt)t>o to Equation (|3.ip . 



To prove Theorem HI we need the following basic result on standard SDEs with jumps. This 
is an easy consequence of Hypothesis Ei For the sake of completeness, we give its proof in 
Appendix iBl 



Proposition 5. For Q G ViE), we set for t G [0,r],x G C 



M 



ip^{t,x) = "^[^l- -1^ 1 when Q{Xt = x) > and ip^{t,x) = f otherwise. 

Let Hypothesis\^hold. Then, for any (xo,yo) £ ^^M x there exists a unique strong solution 
(Xt, )jg[o,T] the following SDE with jumps: 



17 ^«J(T'=-,X^,_) J (3.3) 

yt = yo + f!,b{s,Xs,Ys)ds + f^a{s,Xs,Ys)dWs + fl, -f{s-, X,_,Ys-)dXs. 

Proof of Theorem^ The proof of Theorem S] is split in three main steps. 

• First, we show the existence of a probability measure Q G V{E) solving the Martingale 
Problem (|3.2p . This result is obtained by considering the associated interacting particle 
system: we show that each particle converges in law, and that any probability measure 
in the support of the limiting law solves the Martingale Problem. This is done in 
Section 13.31 



Second, we show the uniqueness of the probability measure Q G V{E) solving ()3.2p . 
To do so, we introduce a function ^' : V{E) — )• V{E) defined as follows. Let {Xq,Yq) 
be a random variable distributed according to Cmit under P. Then, we know from 
Proposition [5] that there exists a unique process {X'j^,Y^)i^[Q^']r] solving 



= Y, + b{s, XQ, Y,^)ds + /* a{s, X^, Y,^)dWs + 7(^-, , 



(3.4) 
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and we define 

^{Q) = law{{X^,Y,%^^,^r])^nE). 

As we will see in Section [3.41 ^ (or more precisely ^ iterated /c-times) is a contraction 
mapping for the variation norm. Combining this result with the following Lemma gives 
the uniqueness of the probability measure solving (j3.2p . 



Lemma 6. Let Q G ViE). We have *&(Q) = Q if and only if Q solves the Martingale 
Problem ()3.2p . In this case, {X'j^,Y^)f^[(j^x] solves Equation ()3.ip . 



Then, the existence and uniqueness of Q satisfying ^(Q) = Q and Lemma [6] au- 
tomatically give the strong existence and uniqueness of {X, Y) satisfying (jS.ip since 
we necessarily have E.[f{Yq^k_)\Xrpk_] = KQ[f{Yj'k_)\Xrpk_] and thus {Xt,Yt)t^[Q^T] = 

[^t ' -^i )te[o,T]- 



Proof of Lemma\^ The direct implication is clear since ^'(Q) = Q gives that E[/(1^^)|XP] = 
ip^{t, X"!^). Thus, (X^, l^'^)(g[Q solves ()3.ip and in particular Q solves the Martingale 



Problem (|3:2|) . 

Conversely, let us assume that Q solves the Martingale Problem ()3.2p . We know from Proposi- 



tionj5]that strong uniqueness holds for the SDE {X^^, ^ )te[o,T]- From lLepeltier and Marchal 
(|l976l . Theorem II13), we know that strong uniqueness implies weak uniqueness, which pre- 
cisely gives ^'(Q) = Q since Q solves the Martingale Problem ()3.2p and therefore the Mar- 
tingale problem associated with ()3.4p . In particular, we have E[/(yj'^)|Xj^] = (p'^{t, X^^), and 
te[o,T] solves (|3.ip . ■ 

3.2 The interacting particle system 

We assume that the probability space (fi, J-", P) carries all the random variables used below. 
Now, we set up the particle system related to the Martingale Problem ()3.2p . In the following, 
will denote the number of particles, (XQ,yg),i S N are independent random variables 
following the law Cmit under P, and {Wf,t > 0),i G N are independent standard Brownian 
motions. We build an interacting particle system {{Xl'^ ,Y^'^),t > 0)i<i<N with the following 
features. For i = 1, . . . ,N, {X^'^ , t > 0) is a Poisson process with intensity: 



A(t-, x-i^)/(y,!i^) Ef=ii 



(3.5) 



and Y^'^ solves the following equation: 



y/'^ = yj + [\{s,Xi''',Y:'^)ds+ fa{s,Xl'^,Y:''')dWi+ f\{s-,r^^X'-)dXi^''- 
Jo Jo Jo 

(3.6) 

In fact, we can give an explicit construction of this particle system, which will be useful 
later and we explain now. Let us consider (?7*''^)j,fceN* a sequence of independent uniform 
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variables on [0, 1] and (-E'*'^)j,fceN* a sequence of independent exponential random variables 
with parameter These variables are independent, and independent of the previously 
defined Brownian motions. We define the times 

k 
1=1 

and we can order {T^''',i = 1, . . . , N,k > 1), such that < T^i'^'i < • • • < T^^^'"' < ... almost 
surely. 

Up to the first jump of X^'^ , ^t'^ defined as the unique strong solution of 

y/'^ = yj+ fh{s,XlY:''')ds+ f a{s,Xi,X''')dWl. 
Jo Jo 

At time r = T*'''^', the process X^''^ makes a jump of size 1 if 

. A(r-, xi'f)f{Y:-'') Ef=i h^.^.^H.-. 

and does not jump otherwise. If a jump occurs, we set Yt'''^ = Y^[i^ + 7(r— , 1"^*^^^) 
and, up to the next jump of X*''^, we define Y^*''^ as the unique strong solution of 

y/"^ = y;''^ + J' b{s, xi"^, y;"^)(is + J' a{s, x;"^, y;"^)dvr;, t > r. 



3.3 Existence of a solution to fl3.2p 

We follow the analysis carried out by Meleard ( 19961 ). pages 69 and 70. We denote by 

^/ vi.N T^i.N\ 

v^t '^t Jte[o,T] 



= jf X^i^i '^cx*'^ Y^^'^) empirical measure given by the particle system. It is a 



random variable taking values in V{E). We denote by vr^ G V{V{E)) the probability law of 
H^. For TT £V{V{E)), we denote by 



JV(E) 



IV(E) 

the mean of vr. Let F : ii^ — t- R be a bounded function which is continuous with respect to 
the Skorokhod topology. It induces an application V{E) — t- M — still denoted by F with an 
abuse of notation — such that 

= / F{z)^i{dz). 



Since tt^ is by definition the probability law of , we have by using Fubini's Theo- 
rem E(F(/i^)) = /^(^^F(^)7r^((i/i) = jj^F{z)I{-K^){dz). On the other hand, we have 

-^(/^^) = ^EiIi^((^i'^>^/'^)te[o,T])- By symmetry, (X*'^, y/'^)tg[o,T] has the same law 
as (^t^'^, y/'^)te[o,T] and we get that: 

E[F((x/'^,y/'^),g[o,T])]= / F(z)/(vr^)(d^). (3.7) 

Je 
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Lemma 7. The sequence {7:^)^ is tight. 



The proof of Lemma [7] is postponed to Appendix [Dl Now, wWe can consider a subsequence 
TT^fe which converges weakly to 7r°° e V{V{E)). Let q £ W , < si <■■■< Sq < s < t, 
cp G C°''(M2,M) and ffi,...,<7, G C;,(R2 
for (j). We set for QeV{E), 



be bounded functions with bounded derivatives 



F 



+ 



A(n,X„)/(y„) 



EQ[f{Yu)\Xu 

b{u,Xu,Yu)dy<p{Xu,Yu) + ^a\u,Xu,Yu)d^cl){Xu,Yu)}dy)jflgi{Xs^,Ysj] (3 



(</.(x„ + 1, y„ + 7(^x, X^, Yu)) - ^{X^,Yu)) 
q 



1=1 



We have to check that Q 1— t- F{Q) is continuous with respect to Q for the weak convergence. 
Since / is continuous by Assumption (jl.ip . we first notice that when Q„ converges weakly 
to Q, EQ^[f{Yt)\Xt = x] converges to EQ[f{Yt)\Xt = x] when Q( Xt = x) > 0, un less for an 
at most countable set of times t depending on Q. Then, following [Meleaxdl (|l996l l. it comes 
out that if s,t, taken outside a countable set depending on 7r°°, F is 7r°°-a.s. 

continuous. In this case, we have: 



fe— >+oo 



F 



2^00 



V{E) 



By definition of , we have: 

TV q 



1=1 



where 



{b{u, Xl;^, Y:'^)dy^iX^^'^, 1^'^) + la'iu, X^^, l^'^)a^<^(Xi^, Y^^^) 



H ^ Af,^ +7(^^,-^« ))-'P(-^n :W )\\du. 

22j = lfiYii' )1^^J^,N^^.^N^ J 

We observe now that [M*'^, M-^'^Jt = for i 7^ j since, by construction these martingales do 
not jump together and {W^,W^)t = 0. Therefore, we get that: 



E[F(^^)2] = 1e {Mr-M^'^fllg,{Xlf,Y,Yf <C/N, 



N 



1=1 



thanks to the boundedness assumption made on functions gi and (p. It comes out that F{Q) = 
0, TT°°{dQ) almost surely. This holds in fact for any function F given by ()3.8p . provided that 
s,t, si, . . . , Sq are taken outside a countable set depending on Tr°°. Since the process {Xt,Yt) is 
cadlag, this is sufficient to show that any measure in the support of tt°° solves the Martingale 
Problem ()3.2|) . In particular, we get the following result. 



Proposition 8. There exists a measure Q G V{E) solving the Martingale Problem ()3.2 
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3.4 Uniqueness of a solution to (13. 2 p 

Let Q G 'P{E) denote a probability measure solving the Martingale Problem (|3.2p . We know 
that such a probability exists thanks to Proposition [8j We want to show that it is indeed 
unique. To do so, we consider another probability Q G T'{E) and study the total variation 
distance Vt(^(Q) - ^(Q)) between Q and Q over E. For t G [0,T], Q,Q € V{E), we denote 



by 
VfO 



and 



l[o,t] 



their restriction to the paths on the time interval [0,t]. We also set 



the total variation distance between 



and 



l[o,t]- 



Lemma 9. Let r := mi{t > 0, X^^ ^ -^?} denote the first time when X'^ and X'^ do not 
jump together. We have 

Vt(^(Q) - ^(Q)) < 2P(t < T). 



Proof. Let us recall that for any signed measure rj on E, the total variation of rj is given by 

V{v) = v^{E) + V^{E), 
where rj = r]~^ — r]~ is the Hahn-Jordan decomposition of rj. Besides, we clearly have 

Iviv) < viv) < V{v), 

where V{r]) = sup{\r]{A)\, A C E measurable}. 
We have for any measurable set A oi E, 

mxf,YSe[o,n^iX?,Y,%,^o,T^) 

> P((X,^, Y,%^[o,T] G A, {X^, Y,%^[o^] A) 

= P((XF, yS^[o,t] G A) - Fiixf, yS^t] G a, (X?, y,^),e[o,T] G A) 

> F{{xf,YS^[o,T] G A)-F{{X^,Y,%^^o,T] G A). 
By taking the supremum over A, we get 

F{{xf,YSmT] / iX^,Y,%^^,^T]) > t^T(^(Q) - ^ 
which gives the claim. 

Up to time r, we have {X'^,Y^) = (Xl^jY^), and we get that 

1 1 



Hr<t} 



ipQ{s,x9) ^Q{s,xf) 



ds 



is a martingale. In particular, we have 
dF{T < t) 



dt 



E 



i|.>,}A(t,x,^)/(y,^) 



1 



1 



< 



^^t,Xl^) ip^t,X^) 

\^%,xf)-^^{t,xf ' 
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Now, let us observe that ^^{t,x) - ip^{t,x) = %[/(^^)i<x.=^}]-%[/(y')i{x.-^}] + 
W^^lQ^Xt = x) - q{Xt = x)]. Thus, we get that \ip'^{t,x) - ^'^{t,x)\ < 



(|EQ[/(yi)l|x,=.}] - EQ[f{Yt)l{x,=.}]\ + miXt = x)- Q{Xt = x)\) , and there- 
fore (we use here that Q is invariant, and is the law of {X^,Y^)t>o) 

E[|y,^(t,xF) - ip^it,xf)\] < Yl l%[/(>^t)l{x.=.}] - EQ[/(yi)l|x.=x}]| +lmXt = x) 

x&Cm 



To check the last inequality, one has to observe that for a simple nonnegative function /(y) = 
E"=i fi^{A,}, where AiCM. are Borel sets, we have \E^[f{Yt)li^Xt=x}\-^Q[f{yt)'^{Xt=x}\\ < 
127=1 fiMXt = x,Yte Ai) - q{Xt = x,Yt e Ai)\ < JVt{Q-Q). By passing to the limit, 
this property holds for any bounded measurable (hence continuous) function /. 
To sum up, we have for t G [0, T], 

— Jl^ r'^ — 17^ — 

Vri^iQ) - ^{Q)) < 2P(t < T) < 4—^ / y,(Q - Q)ds < 4^TVt{Q - Q), 

j_ Jo I 

since Vo(^(Q) — ^'(Q)) = and s i— ?• Vs{Q — Q) is nondecreasing. Let us recall that the set 
of bounded countably additive measures on E endowed with the total variation norm is a 
A 7^ 

Banach space. If i^—p-T < 1, we get by the Banach fixed point theorem that ^ has a unique 
fixed point that is necessarily Q. Otherwise, we get by iterating that: 

1 k 



t2 ,t 



4M^r 



P Jo - kl 



-VriQ-Q). 



< 1 and ^^^^ is a contraction mapping. Thus, ^''-'^^ has a 



When k is large enough, ^ 

unique fixed point that is necessarily Q. 

This concludes the proof of Theorem SI Besides, using the notations of Section 13.31 we get 
that any convergent subsequence of vr^ should converge to ^^((iQ). This gives the weak 
convergence of tt^ towards ^^((iQ). By (|3.7p . we get that any particle converges in law 
towards Q. 

3.5 Convergence speed towards Q 

Now that we have proved that each particle converges to the invariant probability measure, 
we are interested in characterizing the speed of convergence of the interacting particle system 
towards this measure. This question is of practical importance, since one would like to use 
the following approximation 



1 ^ 



i=l 
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and have an estimate of the error involved. 

First, we need to introduce some additional notations. We consider the same particle sys- 



tem {Xl''^ ,Y^'^)t>o as in Section [3.21 constructed with the random variables T*'*^, C/*''^' and 
{W^,t > 0). With these variables, we construct now the processes {Xl,Y^)t>o as the unique 
solution of 



(3.10) 



[Yi = + b{s, XlYi)ds + a{s, XlWdW^ + /o 7(s-, ^i- , Y^-)dXi. 



By construction, the law of {Xl,Y^)i^[Q^'p] is the invariant probability law Q, since ^ 
By using the same argument as in Lemma El we have: 



< 2P((X^^,y/'^),e[o,T] / (X/,F/),e[o,T]) = 2P(ri < T), 



where = inf{t > 0, X} / X/'^}. We also set for i = 2, . . . , TV, = inf{t > O^Xj ^ X^'^}. 

Proposition 10. Let us assume that Cinit is such that ¥{Xq = x) > for any x G Cm- 

Then, there is a constant K > such that 



K 



(r^ < T) < 



Proof. By construction, the processes X^ and X^'^ may become different at the times T^'^ if 



A/ 



ana > T7 — „w,jv — r; 



-, or conversely. 



1 , fc _ 'J 



Thus, l{rl<t}-/o l{rl>s}^ 



have 



A(6.xi)/(y/) 



ds is a martingale and we 



f 

(r^ <T) = 

A/ 



\{tai)f{n) \{t,x^^^)f{Y,''^) 



(it 



Let us observe that on {r^ > t}, we have (X/, Y"/) = (X^ ' ,Y^ ' ). Therefore, 



(r^ < T) < /E 



L{rl>i} 



1 



1 + Ei=2 l{xf^=Xl} 



dt 



.(3.11) 



Now, we study 



iov X £ Cm, and set 
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When qt{x) > 0, we have 

1 1 + Ef=2 ^{Xl'''=x} 



+ 



< 



< 



+ 



1 
1 



fqtix) 



1 / 



We analyze these two terms in a similar manner. We introduce (Xf +\ yj^+^)tg[o,T] another 
copy of (Xj^, y/)jg[Q , which is independent from all other existing processes. We have: 



< 



< 



, N+l 



i=2 



+ 



N 



N 



E ^{Xl=x 



^ N+l 

«*(^) - 1^ E 1 



AT ^ -{Xl=x} 

i=2 



i=2 

TV 



} ^{Xl'''=x} 



+ 



N 



+ ^E^{-^'<o + 



J=2 



1 

iV' 



(3.12) 



since and Xl' may be different only on {r-' < t}. Similarly, we have 



+ iEl{-<0 + ^- (3-13) 



We introduce 



/ TV+l 7V+1 \ 



This is a random function which is independent from {Xl,Y^)t>o- From (|3.11|) . (j3.12p 
and (jS.lSp . we obtain by observing that and have the same law under P: 



<T)< / ¥.[A{X})] + 1 + 



fj \N 



-E 



+ E 



1\ ^T'2<t 



dt. (3.14) 



First, we observe that E[ = ^ ^ gt{a;)>o f^{fy ^ M + 1. Thanks to the independence 

of A{x) and X^, K[A{xf)\Xl = x] = E[A{x)]. On the one hand, we have: 
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E 



< WE 



- / \ 1 /v^JV+l 1 



j=2 ^{Xi=x}. 



y/qt(x){l-qt{x)) 



On the other hand, we have similarly that: 



E 



1 



< 



E 



N+l 



Finally, we obtain that: 



i=2 



iV 



E[A(x)] < 1 + 



/ 



By using the tower property of the conditional expectation, we get: 



E[AiXl)]<(l + ^ 



fj ' 



since we have J2xec V Qt{x) < \/M + 1 by the Caucliy-Schwarz inequality. From ()3.14p 



1 . f f\ f VM + 1 M + 1 f'^ 



TTlr2<t 



dt . 



(3.15) 



So far, we have not used the assumption qo{x) = P(Xo = x) > for any x G Cm- Since the 

A? - _ll'r_ 

jump intensity is bounded by -j-, we necessarily have qt{x) > e L (70(2;), for t G [0, T]. Thus, 

there exists a constant C > (depending on T and £.jnjt) such that < C for t G [0, T], 
X G £j\f . Since and have the same law under F, this gives 



nr' < D < (1 + Z) , ^l±i , cjf P,.' < , (3.6) 

and we easily conclude by Gronwall's lemma. ■ 

Now, we can have an estimate of the accuracy given by the approximation ()3.9p . We assume 
that F : — 7- R is a bounded measurable function. Then, we know by the central limit 
theorem that 

^ E^((^t'^t')*e[o,T]) - EQ[F((Xi,yi),e[o,T])]l AA(0,a2), 



i=l 



where is the variance of F((Xt, )ig[o,r]) under Q. 

Since F is bounded by a constant > 0, we have by Proposition [10] 



■t Jie[0,T]J 



{t^<T}) 



iV" - J] F((xi, y/),,[o,T] )-^Y. ^((^t'"^' 

j=l j=l 

which converges for the L^-norm to when — )■ +00 for a < 1/2. Combining both results, 
we finally get a lower estimate of the convergence rate. 



i=l 
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Corollary 11. Under the assumptions of Proposition [TU. we have for any bounded measurable 
function F : E and any < a < 1/2, 



AT" 



E^[F{{Xt,Yt)t^[o,T]: 



in probability. 



Remark 12. To prove Proposition \l(A we have assumed that P(Xo = x) > for any x € Cm- 
In fact, the same proof would work if we assumed that there existed xq G Cm such that 
P(Xo > xo) = 1 and ¥{Xq = x) > for any x>xq. 

However, in practice, it would have been nice to treat the case ^{Xq = xq) = 1 for some 
Xq € Cm, since we know at the beginning how many firms have already defaulted. Heuristically 



from (|3.15p . we may hope to have for large N that E 



.^2<t < CF{t^ < t) since X'^ 
< M + 1, and and have the same 



and are asymptotically independent, E 
law. This would be enough to conclude. Unfortunately, despite our investigations, we have 
not been able to prove this formally. However, we still observe a convergence speed of l/\/iV 
when Xq = on our numerical experiments (Section^. 



4 Numerical results 



In this section, we illustrate the theoretical results obtained in the previous sections. Let us 
recall that the local intensity (LI) model is a Markov chain with unit jumps occurring with 
the rate A(t— ,x) and that the SLI model is given by equation (jl.3p . First, we highlight that 
the LI and SLI models have the same marginal distributions but different laws as processes. 
Second, we study the convergence of the interacting particle system and obtain numerical 
simulations showing a central limit theorem. 

For our numerical experiments we consider two different models for the process Y described 
below and we assume that the local intensity A and the function / are given by 



A(t, X 



A 1 



X 

M 

/(x) = (xV/)a7 



where the number of defaultable entities M will be taken equal to M = 125 froni now on. 
This choice of A(t, x) corresponds to M independent default times with intensity We also 
assume that there is no default at the beginning, i.e. Xq = 0. 



IS a 



1. In the framework proposed by iLopatin and Misirpashaev (|2008l ). the process Y 
continuous process satisfying a CIR type SDK 

dYt = K{X{t, Xt-) - Yt) + cj^fYtdWt 



(4.1) 

where k, cr > 0. To sample suc h a process, we use the second-order discretization scheme 
for the CIR diffusion given in lAlfonsil km\^ . 



In the framework of Arnsdorf and Halperin ( 20081 ). the process Y is no more continuous 
and may jump when a new default occurs. The process Y solves the following SDK 



dYt = -aYt \o^{Jt)dt + oYtdWt + -fY^-dXt 



(4.2) 
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where a,a > and 7 > 0. Remember that X only has positive jumps. For discretization 
pm'poses, note that between two jump times of X, log(y) solves the following Ornstein 
Uhlenbeck SDE 

dZt = {-aZt + -a^)dt + adWt. 

Even though we could have sampled exactly in this particular case the Gaussian incre- 
ments of Z, we discretize Y using the Euler scheme on Z in our simulation since we 
would make this choice for more general SDEs on Y. 

The LI model can be simulated very easily using a standard Monte-Carlo approach as it is 
sufficient to know how to sample a Poisson process with intensity A; we do not need any 
interacting particle system. 

4.1 Practical implementation 

In this part, we describe our implementation of the particle system to sample from the dis- 
tribution of {Xt,Yt)Q<:t<T- We recall that the process X has no continuous part and makes 
jumps of size 1. The process Y has a continuous part and may jump at the same times as X. 
We consider a regular time grid of [0, T] with step size h = j^: Sk = kh, < k < D. Assume 
we have already discretized Y up to time s^, the discretization of Y at time s^+i is built in 
the following way: 

• If the process X does not jump between time Sk and time Sk+i we use the increment of 
a standard discretization scheme. 

• If the process X jumps at time s with Sk < s < Sk+i, we proceed in three steps: apply 
the previous case between times and s, integrate the jump at time s and finally apply 
the previous case again between time s and Sk+i- 

This scheme ensures all the Y^ are at least discretized on the regular grid {sO)Si,...,sd}- 

Computational complexity. Studying the complexity is of prime importance when 
proposing a numerical algorithm. 

On the one hand, there are D discretization times at which we recalculate the values of 
each y*, which requires 0{DN) operations. On the other hand, the average total number 
of proposed jump dates (ie. the number of steps in loop [5]) is given by the expectation of 
the underlying Poisson process at time T : NT^. The computation cost of the ratio (|4.3p 

is 0{N) and the complexity is 0{ND + N'^T^). Hence, for fixed model parameters, the 

overall complexity of our approach is bounded by 0{DN + N"^). In practice, N is much larger 
than D, and the most computationally demanding part of the algorithm is the numerical 
approximation of the condition expectation involved in the jump intensity. 
The complexity of the interacting particle approach can be well improved if during the algo- 
rithm we keep track of the following two quantities involved in Equation (j4.3p 

N N 

Dl=Y.f^^i-)hxl_=,} and Nl=Y,hxl_=j} j = 0,...,M. (4.4) 

1=1 1=1 
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Algorithm 1 (The naive particle system algorithm with complexity 0{DN + N"^)-) 
t = //Current date. 

ti = for all i = 1 . . . //Last discretization date for the particle i. 
Sample independently according to the initial law. 

■Sfc = // Last date on the regular grid: Sk <t < Sk-\-i- 
while i < T do 

t' = t + £ (^^N^ / /t! is the potential next jump in the whole particle system. 

while t' > Sk + h do 

Sk = Sk + h 
for i = 1 to do 

Update the discretization of from time ti to time Sk- 

ti — Sk 

end for 

t = Sk 
end while 

/ = uniform r.v. with values in {1, . . . , N}. //Index of the particle which may jump. 
Compute the conditional expectation 



1^1=1 ^{Xi=Xi} 



E= _^^=l}^;=^^ . (4.3) 



R = •~X{t' ,X^)f{Y^)E //Compute the acceptance ratio. 

U = uniform r.v. in [0, 1]. 

\i U < R then //We accept the jump. 

Discretize up to time t' . 

Y^ = Y^ + j{t',X^,Y^) 
= X^ + 1 

ti = t' 
end if 
t = t' 
end while 



Since the processes X are unit-jump increasing these quantities clearly vanish for j > 
max/ Xl_ . 

Let t be the last proposed jump time of the particle system. These two vectors can be easily 
updated at time t' which denotes the next possible jump time. If some ticks of the regular 
grid lie in [t,t'), we set t as the last discretization date in this interval and recompute vectors 
{Dl)j and {N^)j using Equation (|4.4p . This happens D times in the algorithm and can be 
done with 0(M + N) operations as explained in Algorithm [2l 

Case 1: If the proposed jump at time t' is not accepted, there is nothing to compute: 

dI = d{ and N^, = . 
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Case 2: Otherwise, let / denote the index of the particle jumping at time t', we use 
■Di,=Dl-f{Y/); Dl^' = Dt' + f(Yj) ifj=X/ 



t ) 

Nl, = Nl - 1; Ni^^ = Ni+^ + 1 if j = Xl, (4.5) 

nI = Nl; d{, = d{ otherwise. 



Algorithm 2 (The improved particle system algorithm with complexity 0{DN).) 
t = 

= for all i = 1 . . . //Last discretization date for the particle i. 
Sample (X*,y*) independently according to the initial law. 
Set L>' = A^' = for < I < M. 

for i = 1 to iV do //Calculate D and N. We can directly set iV° = N , = Nf{YQ 
£)/ _ jy/ _ Q 1 < / < M when the initial law is a Dirac mass at (0, Yq). 

N^' = N^' + 1, D^' = D^' + /(y'O 
end for 

•Sfc = //Last date on the regular grid. 
while t < T do 

t' = t + £ { //t' is the potential next jump in the whole particle system. 



. L 

while t' > Sk + h do 
Sk = Sk + h 
for i = 1 to do 

Update the discretization of from time tj to time Sk- 

ti = Sk 
end for 

t = Sk 

Reinitiahze Z)' = A^' = for < / < M. 
for i = 1 to A^ do //Recalculate D and N. 

N^' = N^' + 1, D^' = D^' + f(Y') 
end for 
end while 

/ = uniform r.v. with values in {1, ... , A^} //Index of the particle which may jump. 

R = ■^X{t' , X^)f(Y^) //Compute the acceptance ratio. 

U = uniform r.v. in [0, 1] 

if U < R then //We accept the jump. 

N^' = N^' - 1, D^' = D^' - f{Yi) 

Discretize Y^ up to time t' . 

Y^ = Y^ + -i{t' ,Y^) 

X^ = X^ + 1 

N^' = N""' + 1, D""' = D""' + f{Y') 

ti = t' 
end if 
t = t' 
end while 
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One should notice that in cases 2 and 3, the updating cost does not depend on the size of 
the vectors {Dl)j and {N^)j. Using these updating formulas, we can improve Algorithm [T] to 
obtain Algorithm [2j With this new algorithm we only have to compute D full approximations 
of the conditional expectation which has a unit a cost of 0{M + N) and the rest of the time 
we use the updating formulas (j4.5p . which happens on average less than NT^ times. Then, 

the overall cost of this new algorithm is 0{D{M + N) +NT^). For fixed model parameters, 

this complexity reduces to 0{DN). This new algorithm has a linear cost with respect to the 
number of particles. Thus, we managed to propose an interacting particle algorithm with 
the same cost as a crude Monte-Carlo method for SDEs since M is in practice fixed and 
much smaller than N . The CPU times of the two algorithms are compared in the following 
examples. 

4.2 Marginal distributions 



I \ 1 \ 1 1 \ 1 1 . 1 . 1 

1 3 5 7 9 11 13 

(a) Default distribution of the SLI model for Y given 
by Equation (g^ with T = 1, Fo = 1, a = 1, = 
0.3,7 = 1,^ = 2.5. 



(b) Default distribution of the SLI model for Y given 
by Equation (|I?T]) with T = l,Yo = I, k. = l,a = 
0.3, A = 2.5. 



(c) Default distribution of the LI model for T 
1, A = 2.5. 



(d) Binomial distribution function with parameters 
A/, p=l-e-^^/^. 



Figure 1: Comparison of the marginal distributions of the LI and SLI models. The simulations 
use N = 50,000 samples, D = 100 and / = 1/3,/ = 3. 
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In Figures [H we draw the probability distribution of X for the LI and SLI models for both 
processes Y considered in this part. The comparison of these graphs highlights how the LI 
and SLI models effectively mimic their marginal distributions; their probability distributions 
look almost the same. Since the LI model corresponds to independent default times with 
intensity the default distribution at time T is actually the Binomial law with parameters 

M and p = 1 — exp 



A rp 



as we can see on Figure 1(d) 





Algorithm [T] 


Algorithm [2] 


Model of Fig. 


1(a) 




239 


1.51 


Model of Fig. 


1(b) 




229 


1.49 



Table 1: Comparison of the CPU times (in seconds) of the two algorithms with the same 
parameters as in Figure [H 

We compare in Table [J the computational times of the two algorithms and the gain obtained 
by the second approach is definitely outstanding. Algorithm [2] massively outperforms Algo- 
rithm [1] by a factor of 150. Of course, this gain will be all the more important as the number 
of particles increases. 

Given the impressive match of the marginal distributions, we would like to numerically inves- 
tigate the difference between their distributions as processes. To do so, we have computed in 
each model the length of the longest interval during which X does not jump defined by 



r = sup{t G [0, T] : 3u G [0, T - t], X, 



u+t- 



Xu} 



(4.6) 



Note that with this definition, t = T when X does not jump on the interval [0,r]. The 
histograms of r in the LI and SLI models are shown in Figure El we can see that, in the SLI, 
the length of the longest interval without jumps for X can be very small with a probability 
much higher than in the LI model (the l.h.s. of the histogram in the SLI model is fatter than 
in the LI model). This impression is reinforced by more quantitative observations. From 
the data used to plot these histograms, we have computed in Table [5] several values of the 
cumulative distribution function of the length of the longest interval without defaults both in 
the LI and SLI models. These quantities differ sufficiently to be numerically convinced that 
these two distributions do not match. 





p(t < r/4) 


P(r < r/8) 


SLI model 


0.1911 (± 0.0033) 


0.0200 (± 0.0012) 


LI model 


0.1645 (± 0.0033) 


0.0113 (± 0.0009) 



Table 2: Some values of the distribution function of the length of the longest interval without 
jumps for T = 2 and A = 2.5. The SLI model is defined as in Figure 2(a) The simulations 
use 50, 000 samples. The values between braces correspond to twice the standard deviation 
of the estimator. 



4.3 Convergence of the interacting particle system 

When introducing the interacting particle system, we emphasized that it was not only a 
theoretical tool but that it was also of practical interest as it satisfies a strong law of large 
numbers. From a numerical point of view, the efficiency of the particle system depends on the 
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(a) SLI model for Y given by Equation H4.2I) with Yo = (b) LI model 

l,a = l,a = 0.3, 7 = 3. 

Figure 2: Histogram of the length of the longest interval without X jumping for T = 2 and 
A = 2.5. These histograms uses 50, 000 samples. 

rate at which every particle converges to the invariant probability (see Section 13. 5|) . In that 
section, we proved that this convergence rate was faster than N" for < a < 1/2 where 
is the number of particles. Now, we want to study this convergence rate in several examples: 
the first example only involves the marginal distribution of the particle system at maturity 
time, whereas the other two examples require the knowledge of the whole distribution and 
not only the marginal ones. 

Number of defaults distribution. First, we start with a simple example. We want to 
study the convergence of the estimator of ¥(Xt = 3) computed on the particle system. We 
ran 5000 independent copies of the interacting particle systems and we computed the value 
of the estimator for each system. In Figure [3l we can see the centered and renormalised 
histogram of the values obtained for the empirical estimators. The histogram can be 
compared to the density of the standard normal distribution plotted as a solid line on 
Figure [3] and they match pretty well. This result suggests that a kind of central limit 
theorem should hold in practice even though we did not manage to prove such a result. 



Asian option on the number of defaults. For our second example, we consider an Asian 
option on the default counting process whose price is given by 



1 



T 



P = E\^-J^ Xudu-K 

This price P will be approximated using the corresponding particle system estimate Pat, 
where N is the number of particles. We are interested in the limiting distribution of P^. 
Because the process X has no continuous part and only makes jumps of size 1 and Xq = 0, 
the pathwise integral can be rewritten 

^ / X^du = Xt-^ Y1 * 

-'^ t<T s.t. Xt-^Xt 
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Figure 3: Centered and renormalised distribution of the estimation of ¥{Xt = 3) using 
the particle system when Y is given by Equation ()4.2p with T = 1,Yq = l,a = l,a = 
0.3,7 = 1, A = 2.5, / = 1/3,7 = 3 and iV = 10000. The histogram was obtained using 5000 
independent particle systems. 



Hence, there is no need to approximate the integral, it can be computed exactly (up to the 
simulation of X). The example requires to sample the joint distribution of Xt and the sum 
of the default times. 

On Figure HI we have plotted the distribution of Pn after renormalizing and centering. As 
before, the solid line is the standard normal density. We can see that the limiting distribution 
looks very much like the Gaussian distribution, but such an histogram does not enable to de- 
termine the rate of convergence to the limiting distribution. Actually, the rate of convergence 
is given by the decrease rate of A/Var(P/v) = (lE(|PAr — Pp))"*^^^. From a practical point of 
view we do not have access to P, so we have approximated it by the empirical mean P of the 
data set used to build the histogram of Figure [H 

We can see on Figure [5] that the rate of decrease of ^JK(\PN'^^^\^ recalls the shape of 
a negative power function. Then, we have decided to compute the linear regression of 
— |logE(|P/v — -Pp) with respect to log(A^) on our simulations of P/v for N varying from 
100 to 10, 000 with a step size of 100, which gives a set of 100 data. The idea of the regression 
is to write 

-^logE(|P;v-P|') = alog(iV) + /3 + e^ 

and to minimize the series X^^v^Af- "^^^ minimum is achieved for a = 0.5014, f3 = —0.1361 
and the empirical variance of the sequence (e^) is equal to 10~^. This computation yields 
that the rate of convergence to the limiting distribution is \/iV. It ensues from this result 
combined with the analysis of the histogram U] that a Central Limit Theorem with rate \/iV 
should hold. 
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Figure 4: Centered and renormalised distribution of the estimator Pjy of the Asian option 
price using the particle system when Y is given by Equation ()4.ip with T = 2,Yq = l,a = 
1,(7 = 0.3, K = 1, f = 1/3, / = 3 and N = 10000. The histogram was obtained using 10, 000 
independent particle systems. 
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Figure 5: Convergence rate of -v/ElP/v — -Pp w.r.t the number of particles A^. 
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Longest interval without jump. In this paragraph, we are interested in the convergence 
rate of the estimator of the length of the longest time interval with no jump. We recall the 
quantity of interest already defined in Equation (|4.6p 

r = sup{t G [0, T] : 3n G [0, T - t], = X.^} 

and we consider its particle system estimator tn ■ To numerically sample from the distribution 
of r, we need to know the joint distribution of the jumping times of X, which is a prime case 
of a pathwise estimator. 

On Figure [6l we can see the centered and renormalized distribution of tn together with the 
standard normal density function plotted as a solid line. Again, the limiting distribution looks 
very much like a Gaussian distribution. Using these 5, 000 independent particle systems, each 
with 10, 000 particles, we can compute an approximation of r, denote f in the following. 
Now, we can run independent simulations of particle systems with a number of particles N 
varying from 100 to 10, 000. We study the rate of decrease of a/EIttv — fp, which according 
to Figure [7] shows a negative power function shape. If we linearly regress log y^E|tjv — Tp 
against log(A^) we find that we can write 

-^logE(|T^-f|2) = alog(iV) + /3 + e7V 

The linear regression yields a = 0.4865, /? = 1.016 and the empirical variance of the sequence 
(e„) is equal to 0.0004. This regression yields that the rate of convergence to the limiting 
distribution is ^/N, which focuses the existence of a Central Limit Theorem with rate ^/N. 




Figure 6: Centered and renormalised distribution of the estimator f of r using the particle 
system when Y is given by Equation ()4.ip with T = 2,Yq = l,a = 1,0" = 0.3, k = 1, 
/ = 1/3, f = 3 and N = 10,000. The histogram was obtained using 10,000 independent 
particle systems. 
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5 Conclusion 

Local intensity models are wide spread for modelling the default counting process in credit 
risk. Recently, more sophisticated models with a stochastic factor involved in the inten- 
sity have been introduced in the literature. These stochastic local intensity models can be 
automatically calibrated to CDO tranche prices by properly choosing the local part of the 
intensity. This particular choice of the local intensity gives rise to a very specific family of 
SLI dynamics for which we have investigated the existence and uniqueness of solutions. This 
theoretical study has been carried out using particle systems, which turned out to be a clever 
tool for the numerical simulation of such dynamics. We have proved that particle Monte-Carlo 
algorithms based on this particle system approach almost surely converge. The theoretical 
study of the convergence rate enabled us to prove that the almost sure convergence took place 
at a rate faster that for any < a < 1/2. Obtaining a Central Limit Theorem type result 
for such particle systems remains an open question, even though we could highlight such a 
behaviour in all our simulations. Last, we have shown that the interacting particle system 
can be sampled with a computational cost in 0{ND), which is the same asymptotic cost as 
a Monte-Carlo algorithm for standard SDEs. 
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A Proof of Proposition [3] 

The scheme of the proof is the following: 

• First, we prove that ^+ : t,x i — > ^{t,x^) is globally Lipschitz in x. Then, {£") admits 
a unique solution on M+. 

• Second, we prove that the solution satisfies Vt > 0, P{t) > and |-P(t)| = 1- 

Step 1: is globally Lipschitz. Let us prove that for all & E x E there exists a 

constant K such that x) — ^+{t,y)\ < K\x — y\. We have 

Bounding this quantity boils down to bound {i^iYj + 2/(i)A(t, i)(742)*) where 



k>0 



Bound for Ai. We easily get from Hypothesis [T] 

E E EE 1^^^.11(4)+ -(?/fe)+i= E E2K^'cii(4)+-(yi.)+i 

<2 sup \nl^\\x-y\. 

Bound for A2. To bound it, we introduce ±^^%^pV(x^)+ and ifc ^S'^"/,^!, {x))+ in 
order to split (^2)^- in three terms. Each of them is bounded in the following way 



/E£o/(0(^^ 



<i|(4)+-(2/i)+|. 



Then, YlieCu Y.j>a f{j)Kt^i){M)] < A(l + 2j)\x - y\. Combining bounds on Ai and A2, 
weget\^+{t,x)-^+{t,y)\<K\x-y\,whereK = 2 sup 1/^1^1 + 2A (l + 2^) . 

Step 2: the solution is positive with norm 1. Let P{t) denote the unique solution of 
{£"). 



28 



First, we prove that \ft £ M+,V(i,i) G Cm x N, > 0. 

Pj{t) satisfies 

We assume that there exists t > such that Pj{t) < 0. We also introduce u := sup{s < t : 
Pj{s) = 0}. Then, we integrate the above equation between u and t. We get 

The l.h.s. is strictly negative whereas the r.h.s. is non negative. Then, Pj{t) > for all 
t G M+. 

Second, we prove Vt G M+, |-P(i)| = 1- Since P is non negative, (|P(t)|)' = 
E.e£,, E,>o(^0'W- Moreover, E.e^,, E,>o(^0'W = 0. Then, \P{t)\ = |P(0)| = 1. 

B Proof of Proposition [5] 

In fact, we can explicitly describe the solution of (j3.3p as follows. We have (Xq, Yq) = (xq, yo). 
Then, up to the first jump of X, Yj is necessarily defined as the unique strong solution of the 
following SDE 

^t = yo+ I b{s,xo,Ys)ds+ [ a{s,xo,Ys)dWs, t<T. 
Jo Jo 

The process {Xt,t > 0) only increases by jumps of size 1, and has at most M jumps since 
X(t, M) = 0. These jumps may only occur at the times T" and a jump do occur at time 
if the following condition is satisfied 



Observe here that for any x G Cm, t i— )• ip'^{t,x) is cadlag since the process {X,Y) is cadlag 
under Q, and (p'^{T"'—,Xt"~) is well defined. If the jump occurs, we have 

Yx" = Yxn— + ^{T^ — , ^T"— )) 

and Yt is defined up to the next jump of X as the unique strong solution of the SDE 

Yt = Yt^ + r h[s,XT„,Y,)ds + r a[s,XT„,Y,)dWs, Tn<t<T. 

J Tn, J Tn 
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C Proof of E[supi^[o,r] l^/'^ll < oo 



We assume that y^'^ satisfies ^M) and X^'^^ is a Poisson process with intensity ()3.5p . Then 



p 




+ 


i 



+ 



<4r'^Uy^o+T'^' / |6(.,Xl'^,y/'^)f + sup 



since X ' jumps at most M times. The r.h.s being increasing w.r.t t, we can replace in the 
l.h.s. |y/' |P by sup„<;j \Yu' 1^. Burkholder-Davis-Gundy inequaUty yields to 



E 



sup [ a{s,Xl'^,Y}'^)dWs 

u<t Jo 



< CpTP/^-'E 



Jo 



ds 



On the other hand, E 



1,N ^rhN-, 



ds 



since the jump intensity of X^'^ is upper bounded by Now since b, a, and 7 have a 
sub linear growth with respect to y (see Hypothesis ([2])), we get 

E[sup <c( r 1 + E[sup |yi'^|P](is) , 

which gives the result by Gronwall's lemma. 



D Proof of Lemma [7] 



The proof of this lemma is done in two steps. First, we claim that {tt )n is tight if and only 
if the sequence {{Xl'^ ,Y^''^)^^^q j'^) ]^ is tight . To c heck this, we first notice that 'P{E) is a 



Polish space. By Proposition 4.6 in Meleard ( 19961 ). (7r^)7v is tight if and only if (/(vr^)) 



is tight. Then Prohorov's Theorem (tightness is equivalent to sequential compactness) gives 
with equation (|3.7p the claim. 

Now, we must show that ((-'^Z'^, y/'^)jg[o,T])Af is tight. We use Aldous' criterion. 

First, we have to check that, for any e > 0, there exists a constant K such that 

P(sup^g[Q2^] l^/'^l + > K) < e. This is trivial since xl'^ is bounded and 

E[supig[Q 2-] < 00 (see Appendix [U|) . Second, we have to check that for any £ > 0, 

> 0, there exist J > and uq such that 

sup sup P(|(Xif ,y,^'^) - {Xlf,Y^f)\ >v)<e, 

n>n(} Ti,r2e7[o_T]:''"i<''"2<Ti+<5 

where 7|o,t] denotes the set of stopping times taking values in [0,T]. We take |(3;,2/)| = 
max(|x|, |y|) and we assume without loss of generality that < 77 < 1. For convenience, we 
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introduce i^^'^ = Yl'^=i '^{T''''<t}- this is a Poisson process with intensity We distinguish 



the two cases: Uri^ > Vr" + 1 (X^'^^ jumps in (ti,T2] ) and v. 
jump in {ti,T2] )■ We have: 

mxlf,Y^^'') - {Xlf,Y^;'')\ > V) <P {i^lf > vlf + 1) 



1,N 



1,N 



IK 



l,N 



1,N 



{X^'^ does not 



i^kn ■■=Pi + P2- 



Since P(z^r2^ > u^l") < E[iy. 



l,N 



^E[r2-ri],wegetPi<5^. 



P2 is bounded by P(/^^ b{s,Xs,Ys)ds + /^^ cj(s, X^, > r/). Using Markov's inequahty, 



we get 



P-? < -^E 



l,N vl,N\ 



Moreover 
E 



b{s,Xl''^,Yy)ds+ I ais,X^'\Y,^^'^)dWs 

2" 



n 



< 2 E 



2~ " 2~ \ 

l\{s,Xy,Yy)ds +E l\{s,Xl'^,Yy)dWs j 



On the one hand, we have 



b{s, Xi'^, Ys^'^)ds < 6 C(l + \Ys'''' for some con- 
stant C > by Hypothesis [2j On the other hand, Burkholder-Davis-Gundy's inequahty gives 



E 



< CE 



Since E[supj<;2^ l^/'^P] < 00, we 



get P2 < Thus, combining the upper bounds on Pi and P2, Aldous' criterion is satisfied 
for 5-^- 
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